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ABSTRACT 



We present a study of analytic models of starless cores whose line profiles have 



'infall asymmetry," or blue-skewed shapes indicative of contracting motions. We 



compare the ability of two types of analytical radiative transfer models to re- 
produce the line profiles and infall speeds of centrally condensed starless cores 



whose infall speeds are spatially constant and range between and 0.2 km s 1 . 
The model line profiles of HCO+ ( J = 1 -> 0) and HCO+ ( J = 3 -> 2) are pro- 
^ ! duced by a self-consistent Monte Carlo radiative transfer code. The analytic 

models assume that the excitation temperature in the front of the cloud is either 
■rj- ! constant ("two-layer" model) or increases inward as a linear function of optical 



depth ("hill" model). Each analytic model is matched to the line profile by rapid 
least-squares fitting. 

The blue-asymmetric line profiles with two peaks, or with a blue shifted peak 
and a red shifted shoulder, can be well fit by one or both of the analytic mod- 
els. For two-peak profiles is best matched by the "hill5" model (a five pa- 
rameter version of the hill model), with an RMS error of 0.01 km s _1 , while the 
"twolayer6" model underestimates the infall speed by a factor of ~ 2. For red- 
shoulder profiles, the HILL5 and TWOLAYER6 fits reproduce infall speeds equally 
well, with an RMS error of 0.04 km s _1 . The fits are most accurate when the 
line has a brightness temperature greater than 3 K. Our most accurate models 
tend to reproduce not only the line profile shape, but also match the excitation 
conditions along the line of sight. A better match to the line profile shape does 
not necessarily imply a better match to the infall speed and provide guidance on 
how to minimize the risk of obtaining a poor infall speed fit. 

A peak signal to noise ratio of at least 30 in the molecular line observations is 
required for performing these analytic radiative transfer fits to the line profiles. 
Moderate amounts of depletion and beam smoothing do not adversely affect the 
accuracy of the infall speeds obtained from these models. 



Subject headings: stars: formation — radio lines: ISM — radiative transfer 
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1. Introduction 

There are four distinct stages in isolated star formation: (1) formation of a gravitation- 
ally bound core withing a molecular cloud, (2) the gravitationally driven collapse of that 
core, (3) the formation of a central star-like object, (4) dispersal of the remaining cloud 
material (Shu, Adams, & Lizano 1987). The initial two stages are predominantly studied 
through millimeter- and submillimeter-wave-length observations. The molecular line rota- 
tional transition emission visible in these wavelength regimes allow one to probe both the 
physical and chemical composition of the molecular cloud cores out of which stars form. 
Spectral analysis of the molecular line emission also allows one to investigate the kinematic 
motions within the molecular cloud, particularly the phase of gravitational collapse (Myers 
et al. 2000) as well as the associated bipolar outflows which are ubiquitous around young 
stars (Bachiller 1996). 

In a dense cloud core, commonly used molecular line tracers of moderate optical depth, 
such as CS (J = 2 — > 1), or HCO + (J = 1 — > 0) tend to become self-absorbed. As the cloud 
contracts under the influence of gravity, the central regions tend to become more dense than 
the outer regions. This sets up a gradient in the excitation temperature of the molecular gas 
within the cloud such that the excitation temperatures in the central regions of the cloud are 
greater than those in the outer regions, even if the kinetic temperature is constant within the 
molecular cloud. If the cloud is spherical and static the spectral line emission from this cloud 
in these tracers will be symmetric about the line of sight velocity and self-absorbed. However, 
if there is radial motion around the molecular cloud core, a velocity gradient develops around 
the center of the core, precisely where the excitation temperature is the greatest. Much of 
the emission from the rear of the cloud is not absorbed as it travels through the core, due to 
the high excitation temperature within the core, and continues unabsorbed through the front 
section of the cloud because the velocity gradient between the rear and front sections of the 
cloud results in a Doppler shift which is significant compared to the broadening by thermal 
motions. Emission from the front of the cloud however is absorbed by nearby molecular gas 
with lower excitation temperature and only a small Doppler shift. In the case of inward radial 
motion the rear emission is blue-shifted and the front emission is red-shifted, yielding the 
"classic" blue-asymmetric infall line profile. This line profile has become the predominant 
tool for investigating infall motions within star-forming molecular cloud cores (Snell & Loren 
1977; Zhou 1992, 1995; Mardones et al. 1997; Gregersen et al. 2000; Lee et al. 2001). 

Although this blue-asymmetry results from infall, it is not unique to that phenomenon. 
The superposition of multiple clouds along the line of sight can also cause an asymmetric 
line profile. By observing an optically thin molecular line, which should remain symmetric 
around the line of sight velocity of the cloud core, it is possible to reduce the likelihood 
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of mistaking separate clouds for an infalling cloud. Asymmetric line profiles can also arise 
in kinematically complex regions. Studies have shown that rotation and outflow can also 
produce asymmetric line profiles (Snell et al. f980; Adelson & Leung 1988; Narayanan et al. 
2002). Even the asymmetric line emission of B335, once thought to be the prototypical 
example of spectrally detected inside-out collapse (Zhou et al. 1993, 1994; Zhou 1995), is 
now also thought to also contain an outflow which gives rise to the line wings previously 
modeled as part of the collapse (Wilner et al. 2000). In spite of the multiple molecular cloud 
configurations that can result in blue-asymmetric line profiles, they are still the best way to 
study infall motions available, provided we understand the context in which they occur and 
make clear the necessary caveats. 

Understanding how blue-asymmetric line profiles arise from infall motion is easy, but 
creating a complete model matching the molecular emission from one or more tracers is a 
difficult and time-consuming task. Typically one constructs an infalling physical model and 
uses LVG, Monte Carlo (Choi et al. 1995) or microturbulent (Zhou 1995) radiative transfer 
to model the emission. By iteratively changing the model parameters one is able to match 
some or all of the observations reasonably well. This process can take several minutes per 
iteration, and becomes prohibitively slow as the number of physical parameters becomes 
large and as the number of molecular species increases. The derived cloud parameters are 
also highly dependent on the input physical model. Changing the physical model can lead 
to an entirely different set of solutions. 

"Starless cores" are useful regions for observations of inward motion because they are 
numerous, nearby, relatively simple in their density structure, and their motions are not con- 
fused by outflows. They have been the subject of several studies (Lee & Myers 1999; Lee et al. 
1999, 2001; Bacmann et al. 2000; Murphy & Myers 2003; Tafalla et al. 2002, 2004; Tafalla 
& Santiago 2004; Park et al. 2004), some of which have used radiative transfer models to fit 
their line profiles. In this paper we present simple analytic radiative transfer models, based 
on the simplest available assumptions of possible excitation temperature trends within star- 
less molecular cloud cores. These models can easily be fit to blue-asymmetric self-absorbed 
molecular line profiles, and produce an estimate of the infall velocity suggested by those line 
profiles. We also test the efficacy of the analytic radiative transfer models by comparing the 
infall velocities derived using those models with actual infall velocities in rigorously mod- 
eled condensing molecular clouds. Our main finding is that our analytic radiative transfer 
models can yield good estimates for the physical conditions of the starless cores with few 
assumptions and little computational investment. 
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2. Analytic Models 

We compare two analytic models in this paper. Although it is possible to construct more 
analytic models using the same techniques in this paper we feel that these two models are 
the most appropriate for modeling infall. It is possible to analytically integrate the equation 
of transfer 

'/,., • .7(7 ). (1) 



where Tb is the brightness temperature, T is the excitation temperature and J(T) = 
(hu/k) [exp {hv / kT) — l] -1 . The brightness temperature is defined as Tb = (c 2 /2u 2 k) l v 
and is directly proportional to the specific intensity l v . The general solution to the equation 
of transfer, assuming that optical depth increases away from the observer is 

T B = T ie - T0 + 



[ T0 J(T)e- T dr, (2) 
Jo 



where T is the excitation temperature of a region which varies over the optical depth interval 
from to To, and Tj is the incident specific intensity of radiation on that region at r in units 
of brightness temperature. This equation can easily be integrated over regions of constant 
excitation temperature. We assume regions of constant excitation temperature along the 
line of sight in our first model. 

The first model, first discussed in Myers et al. (1996), we call the "two-layer" model. 
It is perhaps unfortunate that this name has been given to this model and has stuck, as 
it implies that this is a model of two plane-parallel layers along the line of sight. This is 
not true. The model actually applies to a line of sight in which two regions with differing 
excitation temperatures are moving toward each other. We assume the near region has a 
lower excitation temperature than the far region. We assume both regions have a velocity 
dispersion for the observed molecule a and a total optical depth of r at line center. The 
excitation temperature of the front region is Tf, while that of the rear region is T r , and the 
regions are approaching each other with a speed of 2v in . The optical depth of each region at 
a velocity v is 

T A V ) = T exp[-(v - v LSR - v in ) 2 /2a 2 ] , (3) 
T r (v) = r exp [-(v - v LSR + v iQ ) 2 /2a 2 ] , (4) 

assuming that the average line of sight velocity of both regions is ^lsr- Figure 1 provides 
a graphical representation of the model. The brightness temperature of the spectral line, 
obtained from the equation of radiative transfer, is 

AT B (v) = J{T f ) [I - e- T f (v) ] + J(T r ) [l - e - Tr(v) ] e~ T ^ v) - J{T b ) [l - e-^M-^/M] , (5) 
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where is the background temperature. The parameters r , <r, T/, T r , wlsr, and t> in are 
free parameters which can be adjusted to fit blue-asymmetric line profiles. Lee et al. (2001) 
has used this analytic model to derive infall velocity estimates for 29 starless cores. They 
find that infall velocities are typically on the order of a tenth of a kilometer per second, 
comparable to the velocity dispersions measured in these sources. 

Equation 2 is also analytically integrable if J(T) is a linear function of the optical depth 
r. Assuming a simple linear function J(T) = J x + [(J 2 — Ji) /t ]t we can integrate the 
equation of transfer to obtain 

1 - e~ T() 

T B = T,e- T0 + ( J 2 - Ji) + J! - J^- 7 " . (6) 

In the above equation r is a function of the Doppler velocity, so in order to calculate the in- 
tensity of emission across a line profile J(T) must be a linear function of r at every frequency. 
It is possible to construct a series of regions along the line of sight where J(T) varies linearly 
with r within each region and analytically calculate the brightness of radiation at a given 
frequency emitted from that line of sight using the above equation. We have constructed 
our next model using the above technique and equations which we feel approximates the 
excitation conditions seen along the line of sight in infalling clouds. 

The second model, which we call the "hill" model, is introduced for the first time in this 
paper. This model consists of a core with a peak excitation temperature Tp at the center 
and an excitation temperature of T at the near and far edges of the core. The J(T) drops 
linearly from J(T P ) at the center to J(T ) at edges of the core, forming a hill in the J(T) 
profile. The optical depth of the core is r c , and its infall velocity is v c , while the systematic 
velocity of the system is t>LSR- A schematic representation of this model is shown in figure 1. 
We introduce this new type of model because we believe it may be more analogous to the 
excitation profile we are likely to be observing in starless cores. In figure 2 we compare the 
excitation profile as a function of optical depth at line center for a simulated cloud, based 
on a Monte Carlo radiative transfer model, to the best fit hill and two-layer models of the 
line profile derived from that cloud. The excitation profile of our modeled starless core has 
a significant slope, which the hill model can replicate, but the two-layer model, due to its 
use of constant excitation temperature zones, can not replicate. This improvement provides 
motivation for introducing a new analytic radiative transfer model of infall. 

In order to solve the equation of radiative transfer we separate the hill model profile 
into two regions along the line of sight: (1) The portion of the cloud in which the excitation 
temperature is rising along the line of sight, whose optical depth is r/; And (2) The portion 
of the cloud in which the excitation temperature is falling along the line of sight, whose 
optical depth is r r . The optical depth as a function of line of sight velocity for each of these 
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Fig. 1. — Graphical representation of the analytic radiative transfer models utilized in this 
paper. Each graph indicates the excitation temperature for a line transition as a function 
of optical depth at some frequency. Graph (a) is the two-layer model, which has a constant 
excitation temperature Tf for an optical depth of Tf and a higher excitation temperature 
T r for an optical depth of r r . Each region of constant excitation is traveling at velocity i>i n 
toward the other region. The symbol in the left indicates the position of the observer and 
the arrows indicate the direction of motion relative to the observer of each region. Graph (b) 
is the hill model, which has an optional envelope of constant excitation temperature T Q with 
an optical depth of Tp on the side close to the observer, and optical depth on the side 
opposite the observer (shown in dashed lines). In the core of the hill model the J(T) rises 
linearly from J (To) to J(T p ^) over an optical depth of Tf and then falls again to J (To) over an 
optical depth of Tf. In the Rayleigh- Jeans limit J(T) is equal to the excitation temperature 
and the excitation profile is as depicted, however at higher frequencies the lines rising to T p k 
would have some curvature. 
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Fig. 2. — Line excitation temperature versus line center optical depth for three radiative 
transfer models. The solid line indicates the excitation temperature as a function of optical 
depth for the HCO + ( J = 1 — > 0) transition at line center in a starless cloud core whose 
density follows the expression n(r) = n / [1 + (r/r ) a j, with a of 4, and r of 5 x 10 16 cm. 
The cloud has an velocity dispersion of 0.1 km s _1 as well as a constant infall at all radii of 
0.1 km s -1 . The HCO + ( J = 1 — > 0) line profile generated using the Hogerheijde & van der 
Tak (2000) radiative transfer model is shown in the upper right corner of the figure. The filled 
circles indicate the excitation profile of the best HILL5 fit to the emerging line profile of that 
cloud, while the open squares indicate the excitation profile of the best twolayer6 model 
fit to the emerging line profile of this simulated starless core. The line profiles generated by 
these models are shown in the upper right hand corner of the figure. 
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regions is 

Tf (y) = r c exp [-(v - v LS r - v c ) 2 /2a 2 ] , (7) 
T r (v) = r c exp [-(v - v LS r + v c ) 2 /2<j 2 ] , (8) 

assuming a velocity dispersion of a in the entire cloud. We can now solve the equation of 
transfer by integrating along the line of sight through each of these two regions to derive the 
brightness temperature 

AT B (v) = (J(T P )- J{T ))[{l-e-^)/r f {v)-e-^\l-e-^)/r r {v)] 

+ (J(T ) - J(T b )) [1 - e-T-M-T/M] , (9) 

where TJ, is the background temperature. We believe the excitation temperature gradient in 
this model might more closely approximate the excitation profiles in observed clouds. The six 
free parameters in the hill model are rc, cr, T , Tp, i>lsr, and vq. It is also possible to add an 
envelope of constant excitation temperature moving with its own infall velocity. Assuming 
this envelope has an excitation temperature of T , the optical depth of the envelope and the 
velocity of the envelope adds two more free parameters, creating a maximum of eight free 
parameters possible in the hill models. 

In order to reduce the total number of free parameters we consider several variants of 
each model by holding certain parameters fixed. We considered two variants of the two-layer 
model, which we call TWOLAYER5 and twolayer6, as well as 4 variants of the hill model, 
which we call HILL5, hill6, hill6core, and HILL7. The number is each variant name refers 
to the number of free parameters in that variant. The parameterization of the two-layer 
model variants are shown in table 1, while those of the hill model variants are shown in 
table 2. 



3. Physical Models 

The analytic radiative transfer models do an excellent job of fitting blue-asymmetric line 
profiles (Lee et al. 2001). Their efficacy in determining physical parameters of the molecular 
cloud, such as infall velocity, has not yet been studied. In this paper we study how effectively 
we can recover the physical parameters of modeled clouds using both the two-layer and hill 
analytic radiative transfer models. We then use this information to determine the point at 
which the systematic errors of the model dominate over random observational errors. 

In our first three simulations, we modeled molecular clouds with a density law of the 

form 

n(0 = l + (r/r r (10) 
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Table 1. Variants of the two-layer model 



Variant Name 


TO 


^LSR 


Vin 


a 


Tf 


T 


TWOLAYER5 


Free 


Free 


Free 


Free 


T b 


Free 


TWOLAYER6 


Free 


Free 


Free 


Free 


Free 


Free 



Note. - r is the line center optical depth, vlsr is 
the line of sight velocity of the system, v- m is the infall 
velocity of the system, a is the velocity dispersion of the 
molecular species, Tf is the excitation temperature of the 
front layer, and T r is the excitation temperature of the 
rear layer. 



Table 2. Variants of the hill model 



Varient Name 






^LSR 




v c 


a 


T 


Tpk 


HILL5 





Free 


Free 





Free 


Free 


T b 


Free 


HILL6 


Free 


Free 


Free 


Free 





Free 


T b 


Free 


HILL6CORE 


Free 


Free 


Free 





Free 


Free 


T b 


Free 


HILL7 


Free 


Free 


Free 


Free 





Free 


Free 


Free 



Note. — te is the line center optical depth of the optional envelope, 
tq is the line center optical depth of the cloud, v lsr is the line of sight 
velocity of the system, ve is the infall velocity of the envelope, vc is 
the infall velocity of the cloud, a is the velocity dispersion of the 
molecular species, T is the excitation temperature of the edge of the 
cloud as well as the optional envelope, and T p k is the peak excitation 
temperature in the cloud. 
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where no is the central density, r is the turnover radius, and a is the power law of the 
density profile at large radii (Tafalla et al. 2002). This model reproduces the flat density 
region (at r < r ) and power law region (at r > r ) required to model starless cores. It also 
is in good agreement with the density profile of a Bonnor-Ebert sphere (Tafalla et al. 2004) 
as well as the "Plummer-like" density model of Whitworth & Ward-Thompson (2001). The 
clouds modeled in this paper have peak densities (no) of 6 x 10 4 cm" 3 , a turnover radius 
(r ) of 5 x 10 16 cm, an outer radius of 2 x 10 17 cm, and an a of 4. We chose to model 
HCO + rotational line emission from these clouds, as it is often moderately optically thick in 
starless cores. We assumed that the kinetic temperature of the gas is 10 K throughout the 
cloud. We set the relative abundance of HCO + at 2 x 10 -9 (van Dishoeck et al. 1995) and 
the HCO + /H 13 CO + abundance ratio to 64. Depletion is not included here, however we do 
investigate the effect of depletion in §4.6. There are four velocity laws which we investigate 
in this paper. The first (Simulation A) is constant infall velocity for all radii in the cloud. 
The second (Simulation B) is a constant infall velocity for r > r$ and velocity for r < tq. 
The third (Simulation C) is constant infall velocity for r < r$ and velocity for r > tq. 
These simulations are summarized in Table 3. These velocity laws were chosen because 
they are simple velocity laws to investigate. Our analytic radiative transfer models define a 
characteristic single infall velocity which is analogous to Simulations A, B, and C. We set 
the velocity dispersion to 0.1 km s _1 over the entire cloud and typically allow the velocities 
within the cloud to be up to twice the velocity dispersion. These results should all scale with 
the velocity dispersion and the peak line brightness. In addition to the simulations described 
above, we also simulated a different density profile based on the best fit Bonnor-Ebert model 
of B68 (Alves et al. 2001). The parameters of this simulation (Simulation D) are discussed 
in §4.4. 

We use these simulated starless cloud cores to synthesize HCO + (J = 1 — > 0) and 
HCO + ( J = 3 — > 2) molecular line emission using the Hogerheijde & van der Tak (2000) 
1-D Monte Carlo radiative transfer model. The simulated clouds were divided into 20 spher- 
ical shells, each of uniform density, kinetic temperature, radial velocity, and with a radial size 
of 10 16 cm. The level populations within the first 21 rotational levels of the HCO + molecule 
were calculated using the Einstein-A and HCO + with H 2 collisional rate coefficients from 
Monteiro (1985) and Green (1975). We then integrate along the central line of sight to find 
the J = 1 — > and J = 3 — > 2 molecular line spectral profiles of HCO + because they are fre- 
quently observed, and often have blue-asymmetric line profiles attributed to infall (Gregersen 
et al. 1997; Narayanan & Walker 1998; Narayanan et al. 1998, 2002; De Vries et al. 2002). 
We then fit these synthesized spectra with the two-layer and hill analytic radiative transfer 
models to check the ability of the analytic models to recover physical parameters within the 
modeled cloud. 
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Several algorithms exist for minimizing multidimensional functions. There is always the 
danger of minimizing to a "local" minimum (meaning the lowest value within some finite 
region of the function to be minimized) instead of the "global" minimum (meaning the lowest 
function value over the entire parameter space) and it becomes more difficult to avoid local 
minima as your function goes to higher dimensions. There are several strategies to avoid 
local minima, and we have chosen a hybrid minimization algorithm to do just that. We 
use the Differential Evolution (DE) algorithm of Storn & Price (1997) initially to separate 
the local minima from the global minimum, then we use the Nelder-Mead simplex method 
(Nelder & Mead 1965) to optimize the fit. The DE minimization algorithm is applicable 
to many astronomical problems, and is similar in scope to simulated annealing. DE is 
an evolutionary algorithm in that several parameter sets are randomly selected and these 
are randomly modified over several "generations" until acceptable convergence has been 
achieved. The variability in each parameter is scaled by the variation of the solution sets in 
each generation, so as the algorithm converges the region of the search decreases without the 
need to include a decreasing "temperature" as in simulated annealing. DE is also superior to 
simulated annealing in that it is more likely to achieve convergence, is a simpler algorithm, 
and has fewer adjustable parameters (Storn & Price 1997). 



4. Results 

We made four simulations of infalling clouds and used the Hogerheijde & van der 
Tak (2000) Monte Carlo radiative transfer model to extract HCO + ( J = 1 — > 0) and 
HCO + ( J = 3 — > 2) line profiles from those simulations. In this section we will discuss 
how effective our analytic models are at deriving the input velocities in those simulation, 
quantify the relative contributions of both random and systematic errors in model fits of 
simulated data sets, and apply this model to an observed infall profile in a source which has 
been rigorously modeled using Monte Carlo radiative transfer techniques. 



4.1. Simulation A — Constant Infall 

Simulation A, which has a constant infall velocity for all cloud radii, has asymmetric 
line profile typical of those seen in infalling clouds. The HCO + ( J = 1 — > 0) line profiles are 
shown in figure 3, along with their optically thin H 13 CO + ( J — 1 — > 0) counterparts. The 
spectra in the bottom on the left and right side are for a simulation with no infall. The middle 
spectra on the left and right are of a simulation with an infall velocity of 0.1 km s _1 . The 
top spectra on the left and right are from a simulation with an infall velocity of 0.2 km s -1 . 
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Overlain on the HCO + ( J = 1 — > 0) spectra on the left panel are the TWOLAYER5 and HILL5 
best fit models, while on the right the TWOLAYER6 and HILL6 spectra are shown. Figure 4 
depicts HCO + ( J = 3 — > 2) and H 13 CO + ( J = 3 — > 2) spectra from the same simulations 
along with their twolayer5, hill5, twolayer6, and HILL6 best fit spectral line profiles. 
Two variants are not shown in these fits, the HILL6CORE model best fit in most cases has a 
te which is very close to 0, yielding the same parameters and spectral line shape as the HILL5 
model. The hill7 model always manages a very good fit to the spectral line shape, however 
in some cases, especially in spectra with a small red component that is not well separated 
from the blue component (such as the middle and top spectra in figure 4), the parameters 
of the HILL7 model are not very well constrained. A wide range of infall velocities, optical 
depths, and excitation temperatures in the seven parameter hill7 model can yield the same 
line profile. 

There are several interesting features in these spectral line fits. When the infall velocity 
is greater than the velocity dispersion the H 13 CO + line profiles become double peaked. This 
is not due to self-absorption as the lines are optically thin. It is due to the fact that the 
molecular gas along the central line of sight is traveling at two velocities separated by more 
than the velocity dispersion, resulting in two peaks in the optically thin component. This 
feature is unlikely to be observed in H 13 CO + ( J = 3 — > 2) lines as they are too weak to be 
easily detected. Also interesting is that while the dip between the red and blue peaks of the 
spectral line profile is deep, both the hill variants and the two-layer variants do a good job 
a reproducing the spectral line, but as the red peak drops in intensity relative to the blue 
peak the analytic model fits tend to produce a spectral line with a single blue peak and a 
red shoulder with no local maximum but only a flattening of the slope. The two-layer model 
variants tend to produce this type of spectrum at lower infall speeds than the hill model 
variants (Figure 3). We call the spectra with distinct blue and red peaks "dip" spectra 
and those with only a flattening of the slope, "shoulder" spectra. We also see that none of 
the models can reproduce a spectral line which has a minimal red peak superimposed on a 
broadly asymmetric line (as in the top spectra of figure 4). These models tend to produce 
spectra that are similar in shape to two overlapping gaussians of equal width. If we try to 
fit a spectral line that differs greatly from one that can be obtained by adding two gaussians 
of equal width, then we tend to not get good line fits with any of the analytic models. 

In figure 5 we show the quality of the infall velocity fit attained by the analytic model 
fits to the HCO + ( J = 1 — > 0) line profiles in this simulation (left side) as well as the model 
fits to the HCO + ( J = 3 — > 2) line profiles in this simulation (right side). In this simulation, 
for the entire infall velocity range up to twice the velocity dispersion, the HILL5 model 
tends to reproduce the modeled infall velocity with an RMS systematic error of less than 
0.01 km s- 1 in the HCO+ ( J = 1 -> 0) lines and 0.02 km s" 1 in the HCO+ ( J = 3 -> 2) 
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V (km s 1 ) V (km s 1 ) 

Fig. 3. — Monte Carlo simulated spectra (solid lines) of the HCO + ( J = 1 — > 0) and 
H 13 CO + ( J = 1 — > 0) (multiplied by a factor of 5) emission from Simulation A (Constant 
Infall). The infall velocity in the lowest set of spectra is km s -1 , those in the middle are of 
a cloud infalling at 0.1 km s" 1 , while those at the top are of a cloud infalling at 0.2 km s _1 . 
The analytic model fits are overlaid on each HCO + spectrum. The open squares on the left 
indicate the TWOLAYER5 fits, while the filled circles are the HILL5 fits. On the right the 
open squares indicate the twolayer6 fits, while the filled circles indicate the HILL6 fit. 
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Fig. 4. — Monte Carlo simulated spectra (solid lines) of the HCO + (J = 3 — > 2) and 
H 13 CO + ( J = 3 — > 2) (multiplied by a factor of 40) emission from Simulation A (Constant 
Infall). The infall velocity in the lowest set of spectra is km s _1 , those in the middle are of 
a cloud infalling at 0.1 km s" 1 , while those at the top are of a cloud infalling at 0.2 km s^ 1 . 
The analytic model fits are overlaid on each HCO + spectrum. The open squares on the left 
indicate the twolayer5 fits, while the filled circles are the HILL5 fits. One the right the 
open squares indicate the twolayer6 fits, while the filled circles indicate the HILL6 fit. 
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lines. The twolayer.5 model tends to consistently underestimate the infall velocity by over 
a factor of 2 in both transitions when it produces good line fits. The TWOLAYER6 model fits 
the infall velocity with an RMS error of 0.03 km s -1 in both the HCO + ( J = 1 — > 0) lines and 
the HCO + ( J = 3 — > 2) lines. Generally the twolayer6 model does a better job in cases 
where the spectrum is dominated by the blue component and the relative size or integrated 
intensity of the red component is small. This happens to be the true for cases with higher 
infall velocities in the HCO + ( J = 1 — > 0) and HCO + ( J = 3 — > 2) spectra in this simulation. 
The twolayer6 model also performs well when the excitation temperature of the rear layer 
matches the simulated excitation temperature in the cloud. The HILL6 model does a good 
job at fitting the infall velocity for infall velocities greater than the velocity dispersion, but 
tends to greatly overestimate the infall velocity in cases where it is smaller than the velocity 
dispersion. The hill7 model also tends to reproduce the infall velocity fairly well, but not 
as well as the HILL5 model. 

It should be noted that a good fit to the line profile does not always translate to a good 
fit to the physical parameters, such as the infall velocity, in the simulated clouds. In the 
case of the TWOLAYER6 models, often the dip solution is a good fit to the line profile, but 
a worse fit to the infall velocity than the shoulder solution. The hill7 model is able to fit 
just about any line profile, but often yields worse infall velocities than HILL5 model fits with 
larger minimized y 2 . In summary, the infall velocities for Simulation A (constant infall) are 
most accurately reproduced by the HILL5 model with a RMS deviation of 0.016 km s -1 over 
22 cases. The TWOLAYER6 model is significantly worse with an overall RMS deviation of 
0.030 km s" 1 over 22 cases. 



4.2. Simulation B — Envelope Infall 

Simulation B, which has a constant infall velocity for r > r and no infall for r < r , 
has asymmetric line profiles typical of those seen in infalling clouds. Figure 6 shows the 
HCO + ( J = 1 — > 0) line profiles for three different infall velocities. The self-absorbed solid 
line in the lowest left and right panels is the line profile for a cloud with no infall, as well 
as the optically thin H 13 CO + ( J = 1 — > 0) line profile. In the middle left and right are the 
HCO + and H 13 CO + ( J = 1 — > 0) line profiles for a cloud with an infall velocity of 0.1 km s _1 , 
equal to the velocity dispersion in that cloud. The upper left and right panels are the HCO + 
and H 13 CO + ( J = 1 — > 0) line profiles for a cloud with an infall velocity of 0.2 km s -1 . 
These line profiles are slightly narrower than those seen in Simulation A (Constant Infall). 
This is due to the fact that the cores, where the most dense gas is located, are stationary 
and not producing any broad line wings. As in our previous simulation we have overlaid 
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Fig. 5. — Infall velocity fits. The above figures indicate the infall velocities obtained by 
fitting TWOLAYER5, twolayer6, HILL5, HILL6, and hill7 models to the Monte Carlo sim- 
ulated spectra in Simulation A (Constant Infall). The left hand panels are the results of 
fitting the HCO + (J = 1 — > 0) spectra, and the right hand panels are the results of fitting 
the HCO + (J = 3 — > 2) spectra. The top two figures indicate the infall velocities obtained 
from the two-layer models. Below the two-layer results are the results of the hill models. 
The symbols indicate the number of free parameters, closed circles indicate the 5 free pa- 
rameter models (twolayer.5 and HILL5). Crosses indicate the 6 free parameter models 
(twolayer.6 and HILL6). Open squares indicate the 7 free parameter model (hill7). The 
TWOLAYER6 model often has two local minima, one of which is usually fit by the TWOLAYER5 
model. We show the other local minimum when it is available in this plot. 
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the TWOLAYER5, HILL5, twolayer6, and HILL6 best fit to the HCO+ spectra. Figure 7 
indicates these same results for the HCO + ( J = 3 — > 2) and H 13 CO + ( J = 3 — > 2) spectra in 
Simulation B (Envelope Infall). 

We investigate the efficacy of our models in matching the infall velocity of the Simu- 
lation B (Envelope Infall) clouds in figure 8. The TWOLAYER5 model once again under- 
estimates the infall velocity in all cases and in both transitions, although it monotonically 
increases throughout the simulation. The twolayer6 model tends to yield the same re- 
sult as the twolayer5 model over many of the infall velocities, but jumps to a higher 
infall velocity which in many cases approximates the correct infall velocity (especially in the 
HCO + ( J = 1 — > 0) fits). We investigated these jumps in the twolayer6 model further 
and found that there are two local minima in the y 2 surface for many spectrum shapes. 

The TWOLAYER6 minima we found correspond to "dip" spectra, with two distinct peaks, 
and "shoulder" spectra with one peak and a flattening of the spectrum on the red side. These 
two spectral line fits tend to have very different parameters for the infall velocity which we 
show in figure 9 for the HCO + ( J = 1 — > 0) lines generated in this simulation. The letter with 
the larger size indicates the fit with the lowest x 2 > and the solid line indicates the expected 
value if the infall velocity modeled analytically matches the infall velocity of the cloud. In 
the HCO + (J = 1 — > 0) profiles, the "dip" solution is a better fit to the line profile up to a 
modeled infall velocity of 0.12 km s _1 , however the "shoulder" solution yields a better fit to 
the infall velocity for modeled infall velocities above 0.04 km s -1 . A similar trend is seen 
in the HCO + ( J = 3 — > 2) results. The "shoulder" fits tend to have a higher infall velocity, 
and are often a much better match to the infall velocity. The "dip" spectra are typically the 
equivalent to the twolayer5 model results. The best "dip" spectrum fit of the TWOLAYER6 
model usually results in an excitation temperature of the front layer equal to the background 
temperature and an underestimate of the infall velocity. We therefore recommend that one 
should attempt to find the local minimum corresponding to a "shoulder" fit when using 
the TWOLAYER6 model, and that one can verify that one has found the "shoulder" fit by 
making sure that the excitation temperature of the front layer is significantly higher than 
the background temperature. 

Another trend revealed in figure 8 is that the twolayer6 model does not do a good 
job of fitting high velocity infall using the HCO + (J = 3 — > 2) line profile in this modeled 
cloud. We can see in the upper-most line profiles of figure 7 that the HCO + ( J = 3 — > 2) 
line profiles do not separate into multiple distinct components at high infall velocities in this 
model. As the lines begin to blend into a profile that is not easily separable into multiple 
components, the two-layer model begins to break down and does not accurately match the 
actual infall velocity. 
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Fig. 6. — Monte Carlo simulated spectra (solid lines) of the HCO + ( J = 1 — > 0) and 
H 13 CO + ( J = 1 — > 0) (multiplied by a factor of 5) emission from Simulation B (Envelope 
Infall). The infall velocity in the lowest set of spectra is km s -1 , those in the middle are of 
a cloud infalling at 0.1 km s -1 , while those at the top are of a cloud infalling at 0.2 km s" 1 . 
The analytic model fits are overlaid on each HCO + spectrum. The open squares on the left 
indicate the twolayer5 fits, while the filled circles are the HILL5 fits. One the right the 
open squares indicate the twolayer6 fits, while the filled circles indicate the HILL6 fit. 
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Fig. 7. — Monte Carlo simulated spectra (solid lines) of the HCO + (J = 3 — > 2) and 
H 13 CO + ( J = 3 — > 2) (multiplied by a factor of 40) emission from Simulation B (Envelope 
Infall). The infall velocity in the lowest set of spectra is km s _1 , those in the middle are of 
a cloud infalling at 0.1 km s" 1 , while those at the top are of a cloud infalling at 0.2 km s" 1 . 
The analytic model fits are overlaid on each HCO + spectrum. The open squares on the left 
indicate the twolayer5 fits, while the filled circles are the HILL5 fits. One the right the 
open squares indicate the twolayer6 fits, while the filled circles indicate the HILL6 fit. 
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Fig. 8. — Infall velocity fits. The above figures indicate the infall velocities obtained by 
fitting TWOLAYER5, twolayer6, HILL5, HILL6, and hill7 models to the Monte Carlo sim- 
ulated spectra in Simulation B (Envelope Infall). The left hand panels are the results of 
fitting the HCO + ( J = 1 — > 0) spectra, and the right hand panels are the results of fitting 
the HCO + (J = 3 — > 2) spectra. The top two figures indicate the infall velocities obtained 
from the two-layer models. Below the two-layer results are the results of the hill models. 
The symbols indicate the number of free parameters, closed circles indicate the 5 free pa- 
rameter models (twolayer.5 and HILL5). Crosses indicate the 6 free parameter models 
(twolayer.6 and HILL6). Open squares indicate the 7 free parameter model (hill7). The 
TWOLAYER6 model often has two local minima, one of which is usually fit by the TWOLAYER5 
model. We show the other local minimum when it is available in this plot. 
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Fig. 9. — We find that there are often two local minima in the \ 2 plane of the TWOLAYER6 
fit to a spectrum. We trace each local minimum for the HCO + ( J = 1 — > 0) spectra in 
simulation B here. The "S" points correspond to a local minimum that tends to produce a 
"shoulder" spectrum, while the "D" points corresponds to a local minimum that tends to 
produce a "dip" spectrum. The larger letter indicates which minimum is the global minimum 
for a particular infall velocity. 
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The lower panels of figure 8 show the hill model variants infall velocities compared to 
the infall velocities in those same clouds. Once again we see that HILL5 does a good job of 
matching most of the infall velocities using the HCO + ( J = 1 — > 0) line profile, but tends 
to break down at high velocities. Overall the HILL5 model fits to the HCO + ( J = 1 — > 
match the simulated infall velocities with an RMS error of 0.08 km s _1 over the full range 
of simulated infall velocities, however for infall velocities less than 1.5 times the velocity 
dispersion the RMS error is less than 0.01 km s _1 . The HILL5 infall velocity fits break 
down at an even lower infall velocity in the case of the HCO + ( J = 3 — > 2) line, however 
the overall RMS error in the infall velocity fits for that transition are only 0.07 km s _1 . 
We believe this is because the contribution of the red component is so low in these cases. 
When the magnitude of the red component drops HILL5 does not produce a good profile fit 
(see the topmost spectrum of figure 6). We also believe that the HILL5 model does a poor 
job when the shoulder starts to disappear yielding a sloping asymmetric line such as the 
uppermost spectrum of figure 7. Once again the HILL6 model does a poor job at fitting low 
infall velocities. The HILL7 does a good job at fitting all the infall velocities in both line 
transitions, typically within 0.03 km s~ x of the modeled value. The error in the fit begins to 
increase, especially in the case of the HCO + ( J = 1 — > 0) fit, as the infall velocity increases. 
This is due to the diminishing magnitude of the red-shifted component. As the blue-shifted 
component begins to dominate the hill7 model fit, the seven parameters of the model begin 
to overdetermine the line profile, resulting in a range of solutions that fit equally well. In 
practice it is possible to avoid this by reducing the number of parameters and using the 
HILL6 model, however in most observed infall profiles the magnitudes of both the red and 
blue components are not as large as seen in the top panel of figure 6 (cf. Lee et al. 2001). In 
typical line profiles, such as those seen in this spherically symmetric infalling cloud, with a 
uniform infall velocity in the outer layers, the HILL5 does an excellent job at extracting the 
infall velocity, while the twolayer.6 does an adequate job for lines with a small red peak 
or narrow shoulder component. 



4.3. Simulation C — Core Infall 

We confine infall to the core (r < r ), and surround the infall region with a static 
envelope in Simulation C. The resulting HCO + and H 13 CO + ( J = 1 — > 0) line profiles are 
shown in figure 10. The HCO + and H 13 CO + ( J = 3 — > 2) line profiles are shown in figure 11. 
Once again modeled infall velocity, this time in the core, is increasing upwards. The bottom 
spectra are for no infall velocity, the middle spectra are for an infall velocity equal to the 
velocity dispersion of 0.1 km s -1 , and the top spectra are for an infall velocity of 0.2 km s _1 . 
The first thing we notice is that the HCO + ( J = 1 — > 0) lines are not very blue-asymmetric. 
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Although the width increases as the infall velocity in the core increases, the infall asymmetry 
does not because the J = 1 — > transition is mostly sensitive to motions in the envelope. 
The emission from the core is mostly absorbed by the envelope because the optical depth 
in the J = 1 — > line is high. The HCO + ( J = 3 — > 2) transition, which is more optically 
thin, does a better job of probing the motions in the core, and shows a greater degree of 
asymmetry as the infall velocity increases. 

The low level of asymmetry in the J = 1 — > line results in a relatively poor fit to the 
infall velocity for both the two-layer model variants and the hill model variants. Figure 12 
depicts the efficacy of the model variants at matching the simulated infall velocity in this 
simulation. The twolayer5, twolayer6, hill5, and HILL7 all significantly underestimate 
the infall velocity from the HCO + ( J = 1 — > 0) line profiles. This is not surprising as that 
transition is very thick and not very sensitive to the motions in the core, where the infall is oc- 
curring. The HILL5 fit tends to be consistently a factor of ~ 5 lower than the simulated infall 
velocity. The situation is more encouraging when examining the thinner HCO + ( J = 3 — > 2) 
line profile fits. The twolayer5 model continues to significantly underestimate the infall 
velocity, but the TWOLAYER6 model matches the simulated infall velocity well for velocities 
between 0.1 km s _1 and to 0.15 km s -1 . Overall the RMS error on the infall velocity fits of 
the TWOLAYER6 model to the HCO + ( J = 3 — > 2) lines is 0.03 km s -1 . The HILL5 model 
continues to underestimate the infall velocity from the HCO + ( J = 3 — > 2) line, resulting in 
an RMS error on the fit infall velocity of 0.07 km s -1 . We believe this is due to the fact 
that the red peak quickly becomes a shoulder in this simulation, and the HILL5 model tends 
to do much better when fitting a distinct red peak. The HILL6 model overestimates the 
infall velocity, while the hill7 model tends to do about as well as the TWOLAYER6 model 
at estimating the infall velocity when fitting the HCO + ( J = 3 — > 2) lines in this simulation. 
Both the hill7 model and the twolayer6 model begin to underestimate the infall velocity 
for core infall velocities greater that approximately 0.14 km s -1 . At this point we run into 
the problem that the line profile is not becoming more asymmetric, just wider, due to the 
increasing velocities in the core. 



4.4. Simulation D — Dense Bonnor-Ebert Sphere 

In our fourth simulation we choose a different density profile and a uniform velocity 
throughout the cloud. We specifically choose the Bonnor-Ebert fit to the density profile of 
B68 obtained by using extinction measurements of the background K giant population to 
produce an azimuthally averaged column density profile (Alves, Lada, & Lada 2001). We 
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Fig. 10. — Monte Carlo simulated spectra (solid lines) of the HCO + ( J = 1 — > 0) and 
H 13 CO + ( J = 1 — > 0) (multiplied by a factor of 5) emission from Simulation C (Core In- 
fall). The infall velocity in the lowest set of spectra is km s _1 , those in the middle are of 
a cloud infalling at 0.1 km s -1 , while those at the top are of a cloud infalling at 0.2 km s~ x . 
The analytic model fits are overlaid on each HCO + spectrum. The open squares on the left 
indicate the twolayer5 fits, while the filled circles are the HILL5 fits. One the right the 
open squares indicate the twolayer6 fits, while the filled circles indicate the HILL6 fit. 
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Fig. 11. — Monte Carlo simulated spectra (solid lines) of the HCO + ( J = 3 — > 2) and 
H 13 CO + ( J = 3 — > 2) (multiplied by a factor of 40) emission from Simulation C (Core In- 
f all) . The infall velocity in the lowest set of spectra is km s _1 , those in the middle are of 
a cloud infalling at 0.1 km s -1 , while those at the top are of a cloud infalling at 0.2 km s -1 . 
The analytic model fits are overlaid on each HCO + spectrum. The open squares on the left 
indicate the twolayer5 fits, while the filled circles are the HILL5 fits. One the right the 
open squares indicate the twolayer6 fits, while the filled circles indicate the HILL6 fit. 
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Fig. 12. — Infall velocity fits. The above figures indicate the infall velocities obtained by 
fitting TWOLAYER5, twolayer6, HILL5, HILL6, and hill7 models to the Monte Carlo 
simulated spectra in Simulation C (Core Infall). The left hand panels are the results of 
fitting the HCO + ( J = 1 — > 0) spectra, and the right hand panels are the results of fitting 
the HCO + ( J = 3 — > 2) spectra. The top two figures indicate the infall velocities obtained 
from the two-layer models. Below the two-layer results are the results of the hill models. 
The symbols indicate the number of free parameters, closed circles indicate the 5 free pa- 
rameter models (twolayer.5 and HILL5). Crosses indicate the 6 free parameter models 
(twolayer.6 and HILL6). Open squares indicate the 7 free parameter model (hill7). The 
TWOLAYER6 model often has two local minima, one of which is usually fit by the TWOLAYER5 
model. We show the other local minimum when it is available in this plot. 
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calculate the Bonnor-Ebert density profile using the modified Lane-Emden equation 



where £ = (r / ' a)y/AirGp c is the non-dimensional radial parameter, a is the isothermal sound 
speed, p c is the volume density at the cloud center, and = — ln(p/p c ). We split this 
second order differential equation in to two first order equations by using the variable substi- 
tutions suggested by Ballesteros-Paredes, Klessen, & Vazquez- Semadeni (2003) and integrate 
the equations numerically to derive the density profile. We use £ max = 6.9, the free parameter 
derived by Alves et al. (2001) to be the best fit to the B68 cloud extinction observations. Di 
Francesco et al. (2002) find that the relative abundance of HCO + to molecular hydrogen is 
less than 1.4 x 10 -10 , which is significantly lower than the value of 2 x 10~ 9 which we have 
been using. Under such low abundances the HCO + lines we would derive from radiative 
transfer calculations within this model cloud would be optically thin. Since optically thick 
lines are required for our analysis, we choose to maintain the HCO + to molecular hydrogen 
relative abundance at 2 x 10~ 9 and only change the density profile in this simulation. We 
divide the spherical cloud into thirty equally spaced shells, the central shell has a molecular 
hydrogen density of 2.6 x 10 5 cm -3 , while the outermost shell has a molecular hydrogen 
density of 1.6 x 10 4 cm -3 and a maximum radius of 1.87 x 10 17 cm or 12,500 AU. We assume 
a constant kinetic temperature of 11 K throughout the cloud, as obtained by Lai et al. (2003) 
using NH 3 (1,1) and (2,2) hyperfine line observations conducted at the DSN 34m telescope. 

We assume a constant infall velocity throughout the cloud at all radii. We mod- 
eled 11 infall velocities ranging from to twice the velocity dispersion, which we 
chose to be 0.12 km s _1 . We base our choice of the velocity dispersion on the 
N 2 H+ (J = 1 -> 0) line width in B68 reported by Caselli et al. (2002a). HCO+ ( J = 1 — > 0) 
and H 13 CO + ( J = 1 — > 0) spectra from this simulation are shown in figure 13. As in pre- 
vious simulations we choose a representative sample of three spectra in the direction of the 
cloud center at different infall velocities. In this case the infall velocities are km s _1 , 
0.12 km s^ 1 , and 0.24 km s" 1 . On the left side of figure 13 the best fit twolayer5 and 
HILL5 models are shown, while on the right side the best fit TWOLAYER6 and HILL6 models 
are shown. The first thing to note is that the H 13 CO + lines are slightly asymmetric, similar 
to H 13 CO+ (J = 1 -> 0), HC 18 0+ (J = 1 -> 0), and D 13 CO+ (J = 2 -> 1) observations of 
L1544 seen in figure 1 of Caselli et al. (2002b). This is as a result of self-absorption, because 
although the H 13 CO + molecule has an abundance 64 times lower than that of HCO + , it is 
not so low that it is optically thin. It is important to be aware of that asymmetry may arise 
in less abundant species, however the degree of asymmetry is indeed less than the asymmetry 
of the more abundant HCO + molecule. As in previous simulations, the line profiles gener- 
ated by the analytic models reproduce the Monte Carlo generated profiles very well. Similar 
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results for the HCO+ ( J = 3 -> 2) and H 13 CO+ ( J = 3 -> 2) lines are shown in figure 14. 
The H 13 CO + lines are more symmetric in this case as there is less optical depth along the 
line of sight. 

We compare the infall velocities in this simulation with those obtained by fitting the 
line profiles with our analytic models in figure 15. As in previous sections we present the 
fits to the HCO + ( J = 1 — > 0) line profiles on the left panels of the figure, and those to the 
HCO + (J = 3 — > 2) line profiles on the right panels of the figure. The upper panels refer to 
the two-layer model fits, while the lower panels refer to the hill model fits. Neither of the two- 
layer models reproduce the simulated infall velocity across the entire range of infall velocities 
we tested. The RMS systematic error in the twolayer.6 model fits of HCO + ( J = 1 — > 0) 
is 0.03 km s _1 , and 0.04 km s" 1 with respect to the HCO+ ( J = 3 -> 2) fits. The HILL5 
model does a better job at estimating the infall velocity in the simulated clouds. The RMS 
systematic error in the HILL5 estimates of the infall velocity based on the HCO + ( J = 1 — > 0) 
profile fits is 0.03 km s _1 , while the RMS systematic error based on the HCO + ( J = 3 — > 2) 
profile fits is 0.01 km s _1 . As we have seen in previous simulations the TWOLAYER6 model 
yields a closed fit to the infall velocity at higher infall velocities. The HILL5 model provides 
a fairly good fit for all infall velocities. These results are similar to the results obtained in 
Simulation A (Constant Infall), which is the other simulation in which infall was constant 
at all radii. This is reassuring as the central densities in these two simulations differ by a 
factor of 4 and they have different radial density profiles. This suggests that these analytic 
models may be able to predict infall velocities from optically thick line profiles for a wide 
range of densities and radial density profiles in star forming clouds. 



4.5. Error Analysis 

There are two components to the error in any fit of a model to data. The first, systematic 
error, is error which arise in the model's description of the data. It is due to the fact that the 
model may not correctly describe the nature of the observed region. The second component 
is random error, which arises from noise in the observation. Although random error can be 
reduced over time or with better instrumentation, usually up to some limit, it is always a 
factor in any observation, and in fitting a model to those observations. In the preceding 
sections we examined how well our analytic radiative transfer models fit synthetic noiseless 
spectra. From this we found systematic errors due to the models shown in figures 5 and 8. 
In order to investigate the behavior of the fits in the presence of random error we decided to 
add noise to the spectra. 
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Fig. 13. — Monte Carlo simulated spectra (solid lines) of the HCO + ( J = 1 — > 0) and 
H 13 CO + (J = 1 — > 0) (multiplied by a factor of 2) emission from Simulation D (Dense 
Bonnor-Ebert Sphere). The infall velocity in the lowest set of spectra is km s _1 , those 
in the middle are of a cloud infalling at 0.12 km s -1 , while those at the top are of a cloud 
infalling at 0.24 km s~\ The analytic model fits are overlaid on each HCO + spectrum. The 
open squares on the left indicate the twolayer.5 fits, while the filled circles are the HILL5 
fits. One the right the open squares indicate the TWOLAYER6 fits, while the filled circles 
indicate the HILL6 fit. 
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Fig. 14. — Monte Carlo simulated spectra (solid lines) of the HCO + ( J = 3 — > 2) and 
H 13 CO + (J = 3 — > 2) (multiplied by a factor of 2) emission from Simulation D (Dense 
Bonnor-Ebert Sphere). The infall velocity in the lowest set of spectra is km s _1 , those 
in the middle are of a cloud infalling at 0.12 km s -1 , while those at the top are of a cloud 
infalling at 0.24 km s~\ The analytic model fits are overlaid on each HCO + spectrum. The 
open squares on the left indicate the twolayer.5 fits, while the filled circles are the HILL5 
fits. One the right the open squares indicate the TWOLAYER6 fits, while the filled circles 
indicate the HILL6 fit. 
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Fig. 15. — Infall velocity fits. The above figures indicate the infall velocities obtained by 
fitting TWOLAYER5, twolayer.6, HILL5, HILL6, and hill7 models to the Monte Carlo 
simulated spectra in Simulation D (Dense Bonnor-Ebert Sphere). The left hand panels are 
the results of fitting the HCO + ( J = 1 — > 0) spectra, and the right hand panels are the results 
of fitting the HCO + (J = 3 — > 2) spectra. The top two figures indicate the infall velocities 
obtained from the two-layer models. Below the two-layer results are the results of the hill 
models. The symbols indicate the number of free parameters, closed circles indicate the 5 free 
parameter models (twolayer5 and HILL5). Crosses indicate the 6 free parameter models 
(twolayer.6 and HILL6). Open squares indicate the 7 free parameter model (hill7). The 
TWOLAYER6 model often has two local minima, one of which is usually fit by the TWOLAYER5 
model. We show the global minimum in this plot. 
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Barranco & Goodman (1998) have derived analytic expressions for the contribution of 
random error in a spectrum to the three parameters of a Gaussian fit based upon the work 
of Landman et al. (1982). They find that the random error in the line of sight velocity and 
in the velocity dispersion of a Gaussian fit are the same and equal to 

a, LSR = a CT = 1.06(M V2 ^, (12) 

where ^lsr is the line of sight velocity of the Gaussian, a is the velocity dispersion of the 
Gaussian, 5 V is the width of the spectral channels, T\ is the peak intensity of the Gaussian line 
profile, and <trms is the RMS noise level in the spectrum. Unfortunately deriving expressions 
for the random error in the parameters of the two-layer or hill models is more difficult. We 
therefore use a bootstrap method to investigate the errors in the parameters of these models. 
We generate a noiseless synthetic spectrum from one of the above simulations, add random 
noise to it, and then perform the analytic radiative transfer model fits. We do this 100 times 
with different random noise of the same RMS value in order to record the variance of each 
model parameter. Assuming the errors are normally distributed, we find that the random 
error in each parameter is merely the square root of the variance in that parameter. 

We initially tested the method on Gaussian line fits, whose random errors should follow 
the expression in equation 12. Figure 16 shows our resulting error determinations for different 
RMS noise values (filled circles) versus the expected value (solid line). We express the RMS 
noise in terms of the signal to noise ratio on the lower axis, and also in terms of the integration 
time it would take to achieve a given RMS for a telescope with a system temperature of 200 K 
and a main beam efficiency of 0.5 when observing a Gaussian line with a peak brightness 
temperature of 5 K, along the top axis. We also use a channel width of 0.04 km s -1 in this 
case, less than the input 0.1 km s _1 velocity dispersion of the Gaussian line and typical of 
current spectrometers. The bootstrap estimates are a good fit to equation 12, with minimal 
scatter due to the random nature of the technique. This test indicates that the bootstrap 
method is a robust estimator of the random error in the parameters of the Gaussian fit and 
should be a robust estimator of the random error in the parameters of the analytic radiative 
transfer model fits. 

We apply the bootstrap method to the HCO + ( J = 1 — > 0) line profiles of Simulation A 
(Constant Infall) shown in figure 3. Figure 17 summarizes our findings for the random 
errors in the velocity dispersion, infall velocity, and line of sight velocity parameters of the 
HILL5 model. The closed circles indicate the error in parameters as a function of integration 
time for a spectral resolution of 0.04 km s -1 , while the open circles are those for a spectral 
resolution of 0.16 km s" 1 . The error increases with the lower spectral resolution because the 
self-absorption feature starts to be resolved out at that resolution. The cross in the center 
of each panel indicates the average systematic error obtained from the fits on figure 5. The 
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Fig. 16. — Random errors in a Gaussian fit to the velocity dispersion and line of sight velocity. 
We estimated errors to a Gaussian fit using a bootstrap method and compared those errors to 
the theoretically predicted errors in a Gaussian fit to a noise spectrum derived by Barranco 
& Goodman (1998). The filled circles are our bootstrap estimates for several different RMS 
noise levels. We use signal to noise ratios in the lower x-axis. We also convert the RMS to an 
integration time assuming a velocity resolution of 0.04 km s" 1 in the spectrum, a main beam 
efficiency of 0.5 and a system temperature of 200 K and place that integration time estimate 
on the upper x-axis. The Gaussian we fit in the example has a peak line temperature of 
5.01 K and a velocity dispersion of 0.1 km s _1 . 
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average systematic error in the velocity dispersion is 0.008 km s" 1 , the average systematic 
error in the infall velocity is 0.009 km s" 1 , while the average systematic error in the line of 
sight velocity is 0.005 km s _1 . The solid line in both the velocity dispersion panel and the 
line of sight velocity panel are the Gaussian errors from figure 16. Unsurprisingly the errors 
in the hill model parameters are greater than those in the Gaussian fit, which has fewer 
free parameters of the HILL5 model, however they tend to follow a similar power law trend 
with increasing integration times. What we can see from this plot is that it is possible to be 
dominated by the systematic errors of the hill model after a reasonable amount of integration 
time on a typical millimeter radio telescope. For spectral resolutions adequate to observe 
the line shape, the velocity dispersion random error drops below the systematic level after 
less than 30 seconds of integration, while the line of sight velocity is well determined after 
less than two minutes. Getting to a point where the systematic errors in the infall velocity 
dominate over the random errors takes the longest amount of time, but only about 10 minutes 
for a line whose brightness temperature is approximately 5 K. This means that generally the 
infall velocity is as well determined as it can get by the hill model when the signal to noise 
ratio is greater than ~ 30. We assume, based on our experience with gaussian line fitting 
that maintaining a similar signal to noise ratio in lines of differing peak intensities should 
yield similar errors. This implies an integration time of 30 minutes for a line whose brightness 
temperature is approximately 3 K, and 250 minutes for a line whose brightness temperature 
is approximately 1 K. This means the observations required to obtain the required signal to 
noise can take a prohibitively long time for weak asymmetric line profiles. 



4.6. Depletion and Beam Smoothing 

Depletion in high density regions of a molecular cloud core is a significant effect in 
observed infall regions. Many observations of CO, CS, and HCO + in dense cores show a ring 
distribution thought to result from an abundance drop in the gas phase due to freeze out 
of these molecules onto grains (Kuiper et al. 1996; Willacy et al. 1998; Kramer et al. 1999; 
Alves et al. 1999; Caselli et al. 1999; Bergin et al. 2001; Jessop & Ward-Thompson 2001; 
Bacmann et al. 2002; Tafalla et al. 2002, 2004). We have added depletion in some cases 
to perform an investigation of the efficacy of the analytic models at estimating the infall 
velocity as depletion increases. We adopt the formalism of Tafalla et al. (2002) and set the 
abundance of HCO + to vary exponentially with the density, 

X(r)=X exp[-n(r)/n d }, (13) 

where X(r) is the abundance of HCO + as a function of cloud radius, X is the HCO + 
abundance in the low density limit (2 x 10~ 9 ), n(r) is the density of molecular hydrogen as a 
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Fig. 17. — Systematic and random errors associated with the HILL5 model. We measured the 
average systematic error in the dispersion velocity, infall velocity, and line of sight velocity 
relative to our input models and plot them above as the large crosses. We then added noise 
to simulate integration times between 20 seconds and 600 seconds with velocity resolutions 
of 0.04 km s _1 (filled circles) and 0.16 km s -1 (open circles), and use the bootstrap method 
to measure the random errors in those same parameters. The upper x-axis indicates these 
integration times. The lower x-axis has a measurement of the signal to noise ratio for the 
0.04 km s -1 resolution case. From this it is possible to estimate the signal to noise ratio 
required in order to have the random error drop below the systematic error for a given 
observation. We find that signal to noise ratios greater than 30 are required, using the 
HCO + ( J = 1 — > 0) lines modeled in figure 3 as a template. The solid line indicates the 
theoretical Gaussian errors for dispersion velocity and line of sight velocity as a function of 
time under the same assumptions. 
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function of cloud radius, and is the density at which the abundance drops by a factor of 
e. We included depletion in the Simulation A (Constant Infall) case where v- m = 0.1 km s' 1 
for n d = 10 6 , 10 5 - 5 , 10 5 , 10 4 - 5 , 10 4 , 10 3 - 5 , and 10 3 cm" 3 . We found for the HILL5 model the fit 
infall velocity changed by less that one hundredth of a kilometer per second down to an nd of 
10 4 cm" 3 for the HCO+ ( J = 1 — > 0) line and 10 5 for the HCO+ ( J = 3 -> 2) line. Figure 18 
illustrates that in each of these cases that level of depletion corresponds to a change in the 
line profile from a two-peaked self-absorbed profile to a "shoulder" profile. Severe depletion, 
beyond the levels reported in figure 18 tends to result in a line too weak to be observed. 
Based on this simulation, we infer that moderate levels of depletion which does not render 
the line profile symmetric or unobservable has only a small effect on the HILL5 model's infall 
velocity determinations. 

Another important observational effect is smoothing by a telescope beam. Our simulated 
line profiles were generated by integrating along a pencil beam through the center of a 
simulated cloud. We have generated 2 dimensional grids of pencil beam integrations of our 
simulated clouds and can convolve them to produce a simulated Gaussian telescope beam. 
We begin by projecting the cloud on the sky plane at a distance of 140 pc, so that a beam 
of 60" is equal to a distance of about 0.04 pc in the cloud. We then convolved two cases 
from Simulation A (Constant Infall), where t>; n = 0.1 and 0.2 km s -1 , to generate beam 
sizes of 20" (0.014 pc), 40" (0.027 pc), and 60" (0.041 pc) in the HCO+ ( J = 1 — > 0) and 
H 13 CO + ( J = 1 — > 0) lines. The resulting beam-smoothed line profile seen toward the center 
of each simulated clouds is shown in figure 19. From the left column of figure 19 we see that 
beam smoothing has a similar effect to depleting the core. This is not surprising as both 
depletion and beam smoothing tend to increase the relative contribution of emission from 
the envelope of the molecular cloud. Our investigation indicates that the infall speed is well 
fit in cases where the beam width is less than or equal to the core diameter (10 17 cm in our 
simulation), however this is dependent on both the density structure and velocity structure 
of the molecular cloud. In practical terms, a higher resolution probes a smaller region which 
may have a large infall velocity and therefore a more asymmetric line. 



5. Applying the Analytic Models 

As a final exercise we have applied the HILL5 model to a high signal to noise line profile. 
We have learned from the simulations that the HILL5 model is the most robust model for 
determining the infall velocity of a molecular cloud from the line profile for a majority of 
line profile shapes. In particular the HILL5 model yields a good fit when the line profile has 
distinct blue-shifted and red-shifted peaks with a self-absorption minimum between them. 
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Fig. 18.— The effect of depletion on the HCO+and H 13 CO+ ( J = 1 -> 0) (left) as well 
as the HCO + and H 13 CO + ( J = 3 — > 2) (right) line shapes. The depletion increases as the 
depletion density (rid) drops. In each column we highlight the range of depletions in which 
the transition remains bright (T r >~ 1 K). In most cases the HILL5 model continues to 
provide a good fit to the infall velocity. The only exception is the upper-most set of spectra 
in which the HILL5 model is not applicable according to the prescription given in §5 as the 
dip between the peaks is not sufficiently deep. The H 13 CO + (J = 1 — > 0) spectra have been 
multiplied by a factor of 5, and the H 13 CO + ( J = 3 — > 2) spectra have been multiplied by a 
factor of 20. 
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Fig. 19.— The effect of beam size of the HCO+ ( J = 1 — > 0) and H 13 CO+ ( J = 1 -> 0) line 
profiles for Simulation A (Constant Infall) clouds with infall velocities of 0.1 and 0.2 km s -1 . 
The simulated cloud is placed at a distance of 140 pc from the observer, typical of the nearest 
starless cores, and convolved with various beam sizes typical of current millimeter telescopes. 
The lowest spectra represent a beam size of 20", the middle spectra represent a beam size of 
40", while the uppermost spectra represent a beam size of 60". The left panel is a simulation 
of an infall velocity of 0.1 km s -1 and the right panel is a simulation of an infall velocity of 
0.2 km s -1 . The H 13 CO + ( J = 1 — > 0) spectra have been multiplied by a factor of 5. 
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For line profiles with this shape, in simulations A and D, which had the same infall speed 
at all radii, the RMS error of the HILL5 infall velocity determination is 0.01 km s" 1 , better 
than the overall RMS error in these simulations of 0.02 km s _1 . We recommend using the 
HILL5 model to derive infall velocities from asymmetric molecular line profiles. 

In some cases there is little or no red-shifted peak, but merely a shoulder on the red 
side of the line profile. These cases tend to occur when the infall velocity is greater than the 
velocity dispersion. In Simulation A (Constant Infall) and Simulation C (Core Infall) red- 
shoulder profiles occur in the HCO + ( J = 3 — > 2) transition for simulated infall velocities 
greater than 0.1 km s _1 , in Simulation B (Envelope Infall) red-shoulder profiles occur in 
both the transitions for infall velocities greater than 0.1 km s _1 , while in Simulation D 
(Dense Bonnor-Ebert Sphere) red-shoulder profiles never occur for all the infall velocities 
in our simulations. In these "shoulder" spectra the twolayer.6 model gives a comparable 
infall velocity estimate to the HILL5 model. We note that in Simulation C, where the 
TWOLAYER.6 model performs better than the HILL5 model the red-shoulder model, lines 
tend to be weak (< 1 K). When using the twolayer6 model it is important to choose the 
local minimum (which should also be the global minimum in these cases) that corresponds 
to the "shoulder" solution set. One can verify that the "shoulder" solution has been found 
by examining the best fit excitation temperature of the front layer, which should be greater 
than the background temperature. If the front layer excitation temperature is equal to the 
background temperature then it is likely that the one has found the "dip" solution which is 
usually never a good fit to the infall velocity. 

Why is the HILL5 model good at finding the infall velocity of most clouds, while the 
TWOLAYER6 model is successful only on spectra with an asymmetric shoulder? The key to 
understanding this lies in how well each model matches the excitation conditions along the 
line of sight at each Doppler velocity. In figure 2 we presented cases in which both the HILL5 
and TWOLAYER6 models match the excitation profile of the simulation quite well for r < 3 
at the line center velocity. Although the figure depicts only one Doppler velocity, the match 
continues to hold over the entire line profile. The resulting infall velocities derived from 
these models closely match the infall velocities in the simulated contracting cloud. In cases 
where the infall velocities derived from the analytic models do not match the infall velocities 
in the simulated clouds, the modeled excitation temperature profiles as a function of optical 
depth along the line of sight derived from the spectral line fits do not match the simulated 
excitation temperature profiles. Even in cases where the spectral line match is quite good, 
if the excitation conditions are not well matched then all the physical parameters, including 
the infall velocity, will not match the parameters of the simulated or observed cloud. 

We find that the HILL5 model tends to reproduce the excitation temperature profile 
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as a function of optical depth quite well for a wide range of infall velocities. The reason 
that this is true is that the excitation conditions in a cloud do tend to rise as a function 
of optical depth over the first couple of optical depths, until the excitation temperature 
reaches the kinetic temperature of the core or begins dropping again as the density drops. 
In all of our simulations a linear fit to the rising excitation temperature, starting with an 
excitation temperature equal to the background temperature, is a good approximation of the 
actual excitation temperature profile at all Doppler velocities. As a result the HILL5 model is 
usually sensitive to the physical conditions of the collapsing cloud and yields a good estimate 
of the infall velocity if the cloud radii which comprise the first few optical depths of the line 
profile are infalling. For this reason we recommend using the HILL5 model. Our simulations 
and error analysis suggest that the line profile should have a peak brightness temperature 
greater than 1 K and a signal to noise ratio greater than 30. The line profile must also have 
a well defined dip or shoulder which must be resolved by the spectrometer resolution. If 
these conditions are met the HILL5 model can be expected to provide an estimate of the 
infall velocity with an accuracy of 0.02 km s _1 . Our simulations suggest the accuracy for 
dip profiles alone is 0.01 km s _1 . 

The twolayer6 model tends to reproduce the excitation temperature profile as a 
function of optical depth well only for higher infall velocities which tend to produce spectral 
lines that are very asymmetric or have well defined red shoulders. In these cases it may 
be tempting to use the TWOLAYER6 method to derive the infall velocity, however there are 
several caveats which one must take in to account. In cases with low infall velocities or 
more symmetric line profiles the excitation temperature of the front layer is often low (at or 
near the cosmic microwave background temperature) while that of the rear layers is often 
much higher than the maximum excitation temperature in the cloud, even in cases where 
the TWOLAYER6 model produces a good fit to the emerging line profiles. In these cases 
the infall velocities as other physical conditions derived from the TWOLAYER6 model do 
not match the conditions in the simulated cloud. The reason this occurs is that there is a 
certain amount of degeneracy in reproducing a symmetric or nearly symmetric line profile 
with the TWOLAYER6 model. A very high excitation temperature region behind a very low 
excitation temperature region produces the same nearly symmetric profile as a moderately 
high excitation temperature region behind a moderately low excitation temperature region. 
At higher infall velocities, when the line profile is very asymmetric or has only a shoulder, 
the excitation temperature of the front region can be accurately determined from the red- 
shifted emission of that line profile. This breaks the degeneracy and creates an excitation 
temperature profile along the line of sight that more closely matches the simulation and is 
more well determined overall. As a result the physical parameters of the simulation, including 
the infall velocity, tend to be well determined by the twolayer6 model in these cases. 
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In figure 20 we present a high signal to noise line profiles upon which we will carry- 
out the procedure prescribed above. The line profile is the isolated component of an 
N 2 H + ( J = 1 — > 0) line profile in L1544 observed by Bourke et al. (in preparation). Based 
on our results in the previous sections we choose to model the emission using the HILL5 
model. That best fit is shown by the points in figure 20. We see that the fit is fairly 
good, though it deviates from the line profile by producing a second peak instead of a shoul- 
der. The HILL5 parameters that result in this fit are shown in table 4. In this case our 
HILL5 infall estimate of 0.099 km s" 1 matches well with the best fit TWOLAYER6 estimate 
of 0.095 km s -1 . Unfortunately since we did not create this profile from a model we do not 
know the actual kinematic conditions in the cloud, however there are two other estimates 
obtained using detailed radiative transfer methods. Keto et al. (2004) have fit a model to 
N 2 H + ( J — 1 — > 0) and N 2 H + ( J = 3 — > 2) spectra at 4 positions around L1544. They find 
an infall speed profile that accelerates inwards from the cloud edge to 0.24 km s -1 . Although 
this peak infall speed is much more than our estimate, we note that their fits to L1544 include 
no velocity dispersion. We speculate that the infall speed may be artificially high in order 
to match the line width in the absence of any turbulent velocity component. Bourke et al. 
(in preparation) have observed this particular source in several molecular line transitions, 
include N 2 H+ (J = 1 -> 0), N 2 H+ ( J = 3 -> 2), DCO+ (J = 2 -> 1), DCO+ ( J = 3 -> 2), 
CS (J = 2 -> 1), CS (J = 3 -> 2), CS (J = 5 -> 4), and C 34 S ( J = 2 -> 1). They have used 
the density model of this cloud proposed by Tafalla et al. (2002) and assumed infall at the 
same velocity over the entire cloud along with a kinetic temperature of 10 K to create radia- 
tive transfer models of L1544 and match their observed line profiles. A simultaneous match 
of all the components in the N 2 H + spectral lines currently yields an infall velocity of 0.1 
km s _1 , which is in agreement with our estimate to within a hundredth of a kilometer per 
second. 



6. Conclusion 

We have shown that analytic radiative transfer models are effective at estimating the 
infall velocities of contracting clouds from the blue-asymmetric line profiles of species of 
moderate optical depth in the region of collapse. The models provide good velocity estimates 
with few assumptions. The greatest advantage of the analytic models is that they allow for a 
quick and reasonable estimate of infall velocity with a minimum of computational effort. In 
cases where nothing else is known about the source, they provide the only possible estimate 
of infall, and in regions where a more robust model can be constructed, the analytic models 



Table 3. Density and velocity laws of the simulated clouds 



Simulation Density Law Velocity Law Values 



A 


TMCWC 


v- m = a 






a = 


0.00, 0.02, ... 


, 0.2 km s" 1 


B 


TMCWC 


Vin = | 




a 


if r < r 
if r > r 


a = 


0.00, 0.02, ... 


, 0.2 km s" 1 


C 


TMCWC 


Vin = | 


a 



if r < r 
if r > r 


a = 


0.00, 0.02, ... 


, 0.2 km s- 1 


D 


BE 


v-m = a 






a = 


0.000, 0.024, . 


. . , 0.24 km s- 1 



Note. — The density laws are described in detail in §3. TMCWC refers to the flattened 
power law profile developed in Tafalla et al. (2002). BE refers to the exact solution of 
a hydrostatic isothermal pressure-confined self-gravitating sphere, otherwise known as a 
Bonnor-Ebert sphere. 



Table 4. Parameters of best HILL5 fit to L1544 N 2 H + profile 



To 


^LSR 


^in 


a 






(km s- 1 ) 


(km s _1 ) 


(km s- 1 ) 


(K) 


3.26 


7.183 


0.099 


0.070 


6.14 
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Fig. 20. — Fits of the HILL5 model to a high signal to noise N 2 H + observation of L1544. The 
spectrum shown is the isolated N 2 H + ( J = 1 — > 0) component observed by Bourke et al. (in 
preparation) toward L1544. The filled circles indicates the best HILL5 fit. 
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provide an excellent starting point for further modeling 1 . The main conclusions of this work 
are as follows: 

1. Analytic radiative transfer models can reproduce the blue-asymmetric line profiles seen 
in infalling starless cores. In addition to the two-layer model first presented in Myers 
et al. (1996) we have found a second class of radiative transfer models (the hill models) 
and explored 5 variants of these two classes of analytic radiative transfer models. All 
the model variations were able to produce blue-asymmetric line profiles similar to those 
we observe in starless cores. 

2. Successful analytic radiative transfer models reproduce not only the line profile shape, 
but also match the initial excitation conditions along the line of sight up to or beyond 
an optical depth of 1. When this requirement is met, the modeled infall velocity 
successfully fits the simulated infall velocity. Those models which match the line profile 
shape, but do not match the line of sight excitation conditions cannot reveal the infall 
velocity and velocity dispersion in the cloud. 

3. The HILL5 analytic radiative transfer model is the preferred model when fitting a blue- 
asymmetric line profile. In simulations with constant velocity as at all cloud radii the 
infall velocities attained by fitting the HILL5 analytic radiative transfer model to the 
simulated line profiles match the infall speed with an RMS error 0.02 km s -1 . The 
HILL5 model tends to perform especially well when there are two distinct peaks such 
that the intensity difference between the self- absorption trough and the redshifted peak 
is greater than 10% of the intensity of the blueshifted peak. The RMS error in the 
infall determination of HILL5 fits to simulated spectral lines with those characteristics 
from Simulations A and D is 0.01 km s _1 . 

4. The physical parameters of simulated clouds that produce blue-asymmetric line profiles 
with no distinct shoulder or redshifted peak are generally not well determined by any 
of the model variations we have studied in this paper. We recommend not using 
these models to estimate infall velocities from line profiles that do not have a separate 
redshifted peak or a distinct redshifted shoulder. 

5. A better match to the line profile shape does not necessarily imply a better match 
to the excitation profile or the physical parameters, such as the infall velocity, of the 
observed cloud. This can be seen most clearly in figure 9 where local minimum of the 



1 Implementations of the models used for this paper are available from the authors at http://cfa- 
www.harvard.edu/~cdevries/analyticinfall.html 
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"dip" solution often provides a lower x 2 fit to the line profile shape, but a poorer fit to 
the infall velocity than the "shoulder" solution. Following the instructions in §5 when 
fitting a blue-asymmetric line profile will minimize the risk of obtaining an unusually 
poor estimate of the physical parameters from the analytic radiative transfer model fit 
to the line profile. 

6. All the analytic models lose sensitivity to the infall velocity if there is a significantly 
optically thick static envelope between the observer and the infall region. In these 
cases a more optically thin line can provide a better estimate to the infall velocity 
in the obscured region. Conversely an optically thin tracer may underestimate the 
infall velocity of an infalling envelope around a static cloud core. These methods are 
sensitive to the velocities and conditions in the regions in which the line is becoming 
optically thick. 

7. A peak signal to noise ratio greater than or equal to 30 is required for the systematic 
error on the infall velocity estimate to be greater than the random error. We recom- 
mend achieving a peak signal to noise ratio of 30 in order to minimize the random 
error of a line profile fit using our analytic radiative transfer models, however nothing 
is gained by achieving a higher signal to noise level. Achieving this signal to noise level 
can take a prohibitively long time for lines whose peak antenna temperature is less 
than or equal to 1 K. 

8. Depletion and beam smoothing can greatly affect the observed line profile in infalling 
starless cores. We conclude that, even in the presence of these effects, if the line profile 
is asymmetric and a good candidate to be fit using these analytic radiative transfer 
models (as described in §5), then generally the radiative transfer models do yield an 
accurate estimate of the infall velocity. If the cloud is depleted to such an extent or 
viewed with such a large beam that the line profile becomes very weak or symmetric, 
then these models will not accurately probe the infall region of the starless core. 
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